#include "header.h"

/* -------------------------------------------------------------------------- */
/* Computes all the split flux derivatives. */
void flux_derv(int point,struct node* p,double** var_old,double* divf){
/* -------------------------------------------------------------------------- */
int i,j,neigh;
double dx,dy,dz,theta,phi,**R;
double dvar,si,sj,DeltaPlus[5],DeltaMinus[5];
double var_tilde[5],var_tildeo[5], divp[5];
double fi[5],gi[5],fluxp[5],fluxm[5],flux_yNp[5],flux_yNm[5];
double flux_zNp[5], flux_zNm[5];
double fluxF[5],fluxG[5],fluxH[5],xfluxo[5],yfluxo[5],zfluxo[5];

R = (double**)malloc(3*sizeof(double*));
for(i=0;i<3;i++) R[i] = (double*)malloc(3*sizeof(double));

for(i=0;i<5;i++) divf[i] = 0.0;

xyzflux(var_old[point],xfluxo,yfluxo,zfluxo);
for(i=0;i<p[point].total_conn;i++){ /* Finding tilde variables in the whole stencil */

	neigh = p[point].Conn[i];

	dx = p[neigh].x - p[point].x;
	dy = p[neigh].y - p[point].y;
	dz = p[neigh].z - p[point].z;

	rotationmat(dx, dy, dz, R);

	for(j=0;j<5;j++){
		var_tilde[j]       = var_old[neigh][j];
		var_tildeo[j]       = var_old[point][j];
		DeltaMinus[j] = ( dx*p[point].qx[j]  + dy*p[point].qy[j] + dz*p[point].qz[j] );
		DeltaPlus[j]  = ( dx*p[neigh].qx[j]  + dy*p[neigh].qy[j] + dz*p[neigh].qz[j] );
		}

	for(j=0;j<5;j++){
		dvar       = var_tilde[j] - var_tildeo[j];
		DeltaMinus[j] = 2.0*DeltaMinus[j] - dvar;
		DeltaPlus[j]  = 2.0*DeltaPlus[j]  - dvar;
		var_tilde[j]   = var_tilde[j] - 0.25*( (1.0 - KKK)*DeltaPlus[j] + (1.0 + KKK)*dvar) ; 
		var_tildeo[j]  = var_tildeo[j] + 0.25*( (1.0 - KKK)*DeltaMinus[j] + (1.0 + KKK)*dvar);
		//var_tilde[j]   = var_tilde[j]  - 0.5*( DeltaPlus[j]  ) ; 
		//var_tildeo[j]  = var_tildeo[j] + 0.5*( DeltaMinus[j] ) ; 
		}
	rotate(var_tilde,R);
	rotate(var_tildeo,R);

	/* Computation of the flux divergence starts here */

	yfluxN_plus(var_tildeo, flux_yNp);
	yfluxN_minus(var_tilde, flux_yNm);

	zfluxN_plus(var_tildeo, flux_zNp);
	zfluxN_minus(var_tilde, flux_zNm);

	FluxKinetic(var_tildeo,var_tilde,fluxF);
	//FluxRoe(var_tildeo,var_tilde,angle,fluxF);
	//FluxAUSMDV(var_tildeo,var_tilde,angle,fluxF);

	for(j=0;j<5;j++){
		fluxG[j] = flux_yNp[j] + flux_yNm[j];
		fluxH[j] = flux_zNp[j] + flux_zNm[j];
		}

	for(j=0; j<5; j++){
		divp[j] =  p[point].X[i]*( R[0][0]*fluxF[j] + R[1][0]*fluxG[j] + R[2][0]*fluxH[j] )
		         + p[point].Y[i]*( R[0][1]*fluxF[j] + R[1][1]*fluxG[j] + R[2][1]*fluxH[j] )
		         + p[point].Z[i]*( R[0][2]*fluxF[j] + R[1][2]*fluxG[j] + R[2][2]*fluxH[j] );
		}

	rotate_back(divp, R);

	for(j=0;j<5;j++)
		divf[j] += divp[j] - p[point].X[i]*xfluxo[j] - p[point].Y[i]*yfluxo[j] 
		                   - p[point].Z[i]*zfluxo[j];
	}

for(i=0;i<5;i++) divf[i] = 2.0*divf[i];

for(i=0; i<3; i++) free(R[i]);
free(R);

}
/* -------------------------------------------------------------------------- */
